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p^- Abstract 



As a generic model system of an asymmetric binary fluid mixture, hexadecane dissolved in carbon 
dioxide is considered, using a coarse-grained bead-spring model for the short polymer, and a simple 
spherical particle with Lennard-Jones interactions for the carbon dioxide molecules. In previous 
work, it has been shown that this model reproduces the real phase diagram reasonable well, and 
also the initial stages of spinodal decomposition in the bulk following a sudden expansion of the 



'^ ' system could be studied. Using the parallelized simulation package ESPResSo on a multiprocessor 

O , supercomputer, phase separation of thin fluid films confined between parallel walls that are repulsive 

for both types of molecules are simulated in a rather large system (1356 x 1356 x 67.8 A^, 
^ ' corresponding to about 3.2 million atoms). Following the sudden system expansion, a complicated 

in the early stages the hexadecane molecules accumulate mostly in the center of the slit pore, but 

^^ ' as the coarsening of the structure in the parallel direction proceeds, the inhomogeneity in the 

^^ ■ perpendicular direction gets much reduced. Studying then the structure factors and correlation 

-j-2 . functions at fixed distances from the wall, the densities are essentially not conserved at these 

H , distances, and hence the behavior differs strongly from spinodal decomposition in the bulk. Some 

of the characteristic lengths show a nonmonotonic variation with time, and simple coarsening 

described by power-law growth is only observed if the domain sizes are much larger than the film 

thickness. 



I. INTRODUCTION 



Fluids confined in pores with linear dimensions on the fim to nm scale find increasing 



applications anc 
5|] and dynamic 



are the sub. 
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ect of many studies, both with respect to their static 



iii,y,y,i4, 



111 ] properties. Considering binary fiuid mixtures, it is natural 



to expect that the (enthalpic and entropic) interactions between the pore walls and the fiuid 
particles may differ for both constituents, and then both density and composition develop 
an interesting inhomogeneity in the direction perpendicular to the pore walls. Of course. 



already in the bulk the binary fluid may undergo both vapor- 
phase separation, resulting in complex phase behavior 12, 



iquid unmixing and fluid-fluid 



131]. In thin slit pores, phase 



separation as a thermodynamic phase transition is still possible in the lateral directions 
parallel to the walls jj], and due to the possible interplay with wetting phenomena [l4, la, 
16|,ll7|| complicated phase diagrams are expected even for strictly symmetric mixtures 4l. Il8|. 
Particularly interesting, however, is the kinetics of these phase transitions as a function of 
time t after a quench. For a strictly symmetric binary Lennard- Jones mixture, where one 
species is strongly attracted by the walls, it has recently been shown by Molecular Dynamics 
simulations that the lateral phase separation kinetics is characterized by a power-law for the 



size of the growing domains 
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i{t) oc r, with [29|, 1^ 



a ^ 2/3 if £{t) is in the range of a few 



jennard- Jones diameters. While for simple diffusive 
systems a = 1/3 both in the bulk [19| and in thin fllms, at late enough times 3l|, the 
much faster domain growth seen by Das et al. [29|, |30t| for a conflned fluid binary mixture 
may be due to some hydrodynamic mechanisms, but is not in accord with the theoretical 



expectations [22, 



23 



24 



25 



26 



m 



281]. Thus, it is interesting to study the kinetics of phase 



separation for other models of conflned binary mixtures, in order to clarify which features 
are universal and which features are model speciflc. 

In the present work, we contribute to this problem by studying phase separation for a 
model of a mixture of hexadecane (C16II34) and carbon dioxide (CO2). There are several 
reasons for this particular choice: flrst of all, supercritical CO2 is a very important fluid 



m 



32 



:he chemical industry, useful as a solvent in which various reactions can be carried out 
particularly applications involving polymers. Thus, the system C16H34 -|- CO2 is a 



prototypical polyne. + solvent sy^tenj^. Secondly, a .athe. ..nple coa.e-g.ained .nodel 

for this system has been developed [35| which describes the experimental phase diagram 



rather accurately. Thirdly, spinodal decomposition in the bulk has already been investigated 
for this model by extensive simulations 36|] . It was found that the system is compatible with 
a growth according to i{t) oc t^^^, when i{t) starts to exceed the Lennard- Jones diameters, 
while at late times a crossover to somewhat faster growth occurs. Limitations due to the 
finite linear dimensions of the simulation box preclude strong statements on the growth law 
during the late stages, however. 

In Sec. 2, we shall introduce our model and briefly discuss the simulation technique. In 
Sec. 3, simulation results will be presented and discussed in the light of various theoretical 
considerations. Sec. 4 contains a summary and outlook to future work. 



II. MODEL AND SIMULATION DETAILS 



A. A coarse-grained model for hexadecane + carbon dioxide mixtures 



Although hexadecane is a very short polymer only, an all-atom simulation of hexadecane 
melts would be difficult, since for a simulation study of phase separation kinetics, length 
scales far beyond the size of a molecule need to be explored, and also large time scales are 



mandatory 22 
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281]. Therefore, it is advantageous to use a coarse-grained model. 



Coarse-graining of polymers is usually done by taking a few chemical monomers (CH2 or CH3 
at chain ends, in this case) together into effective monomers, ignoring completely torsional 



potentials [37 



et al. 



38, 



39 



40| . A successful model of this type for C16H34 was proposed by Virnau 



351], incorporating three successive C-C bonds along a chain (plus the corresponding 



hydrogen atoms) into one effective bead, so that a chain of 5 effective monomers is created. 
Effective monomers along a chain are bound together via FENE (finitely extensible nonlinear 
elastic) potentials |4l[ ] 



f/FENE(r) = -33.75£ppln[l - {r/R, 



VPJ 



R 



pp 



1.5(7 



pp 



(1) 



where e^p, cTpp are parameters of the Lennard- Jones (LJ) potential, that acts between all 
beads of the polymer chains (bonded as well as non-bonded ones) 



U{r) = UUr) - f/Lj(rcut) 



U 



LJ 



4e„, 



O-a/^V^ /0-Q/3\6 



(2) 



where a cutoff rcut = 2rmin, rmin = 2^/^aai3 is used and the potential is shifted to zero at 
f = '"cut so that U{r-) is everywhere continuous, with U{r > rcut) = 0. Here a,/3 = p (if 



the particle is an effective monomer of the chains) or a, f3 = s (if the particle is a solvent 
molecule). The parameters £pp, CTpp and ^ss, o'ss are chosen such that the model reproduces the 
experimentally known 42| critical temperatures T^, and critical densities pc for pure C16H34 
and pure CO2, respectively 3J, |35|. Thus, using [42| T^ = 723K and p^ = 0.219g/cm^ has 



yielded SJ, l35| (henceforth we omit the index p) 



5.79 ■ lO^^V, a = 4.52 ■ 10 



-10 



m , 



(3) 



while the experimental results for CO2, T^ = 304K and po = 0.464g/cm^ have yielded |34j.l35| 



0.726e, as 



0.816(7 



(4) 



With these parameters (Eq. ([3]) and (j4])) the coexistence curves in the temperature-density 
plane and the vapor pressures at coexistence as well as the interfacial tension between the 
coexisting phases are reproduced in reasonable agreement with experiment 34, l35|]. An 
even better description of CO2 could be obtained by including the quadrupole-quadrupole 
interaction 43|], but this is out of consideration in the present context. 

The parameters ^ps, cTps for the interactions between CO2 molecules and effective 
monomers are described 3J, |35|] using a modified Lorentz-Berthelot mixing rule 4^ 



o"sp = (o-ss + 0'pp)/2 



'Sp 



S Y ^ss^ 



PP 



(5) 



with 



34 



35| C, = 0.886. While the standard Lorentz-Berthelot mixing rule (^ = 1) would 



yield a phase diagram topology in disagreement with the available experiments 45], Eq. ([5]) 
gives a phase diagram in rough agreement with these experiments 3J, |35|]. In the following, 
we shall choose e = 1 as unit of temperature (taking Boltzmann's constant ks = 1) and 
0" = 1 as unit of length. Fig. [1] shows an isothermal slice through the phase diagram 
at reduced temperature T* = ksT/e = 1.16 in the plane of pressure p and molar fraction 
X = N'^/{N''^ + N^/5) of carbon dioxide, where A^'^ is the number of carbon dioxide molecules 
and A^^ the number of effective monomers of hexadecane. As will be described below, we 
shall simulate pressure-jump experiments where the system suddenly is brought from a state 
in the one-phase region (the initial state is equilibrated at a density pl^^ = pa^ = 0.8 in the 
middle of the slit pore, which would correspond to a reduced pressure p* = pa^/e = 0.34 in 
the bulk system) into the two-phase region by an isotropic increase of the volume available 
to the particles. 



For a system in a thin film geometry, it is also necessary to specify the boundary conditions 
created by the planar walls confining the thin film. We choose an atomistic description of 
these walls, putting particles on a regular (and rigid) triangular lattice of lattice spacing 
o" = 1, in the {x,y) plane at z = 0.01 and z = 0.99Lz, z being the coordinate in the 
direction perpendicular to the walls. The interactions between the wall particles (w) and 
solvent particles or effective monomers are described by the purely repulsive part of the 



LJ-potential, Eq. ([2]), using r, 



cut 



and 



-wp 



1, 



0"„ 



a 



wp 



(6) 



This choice was made to avoid the formation of precursors of wetting layers of one of the 
species, unlike 



29 



30|, where an attractive interaction between the walls and one of the 



species in the binary (A,B) mixture was chosen. We deliberately choose the wall-particle 
interactions symmetric in the present case, to avoid any strong preference of the wall for 

nn 

one of the components in our case. However, since (unlike to [29|, 130|) the present model 
is not a symmetric mixture in the bulk, we do expect some wall-induced concentration 
inhomogeneities for the present model as well. As demonstrated recently for the case of 
colloid-polymer mixtures 46|, an effective attractive interaction due to repulsive walls may 
arise due to purely entropic origin. 



B. Simulation method and preparation of the initial state 



We study the kinetics of phase separation in thin films of CO2 + C16H34 mixtures by 



Molecular Dynamics (MD) methods 47, 



48 



49] . As is well-known, in simple fiuids and 



binary fiuid mixtures hydro d yna mic interactions are important both for the dynamics of 



fiuctuations near equilibrium 
spinodal decomposition 
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50 [and for the kinetics of coarsening in the late stages of 



M, 



25 



26 
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281]. MD simulations (in the microcanon- 



ical NVE ensemble where energy E is conserved for fixed number of particles A^ and fixed 
volume V) include these effects of hydrodynamics implicitly and fully 47|, |48!, |49!]. In fact, 
previous studies of phase separation in the bulk have used this method successfully for both 
simple liquid-yapqr phase separation and for studies of unmixing of binary fluid mixtures 



[see e.g. 



51 



52 



53 



54 



55] ). We apply for our system the software package ESPResSo, 



version 1.9.7h, 2005 56| which is particularly suitable for simulation of coarse-grained soft 



matter systems on parallel computers. 



For the integration of Newton's equation of motion the Velocity Verlet Algorithm 47 



48|, 



49| is applied, choosing an integration time step 6t = 0.002 r, where the MD time unit here 
defined as r = cr(m/£:)^/^ = 1 corresponds to 500 integration steps. The masses m of CO2 
molecules, effective monomers and wall particles are set for simplicity equal to each other, 
and time units are chosen such that m = 1. 

The initial state is created by first using a small simulation box L^ x Ly x L^ with 
L^ = Ly = 20, Lz = 12 (measured in units of a through this paper) with two repulsive walls 
at z = 0.01 and z = 0.99Lz, as described in Sec. 2.1, and periodic boundary conditions in x 
and y-directions. Into this box, hexadecane molecules were inserted, having selfavoiding walk 
configurations, and the CO2 particles were inserted at random position, at molar fraction 
of CO2 X = 0.6, such that the initial state reaches a reduced total monomer density of 
Ptot ~ 0-8 i'^ ^^^ center of the thin film. The CO2 particles were only allowed to be put 
at positions outside of a sphere of radius a = 1 of each bead, to avoid that in the initial 



state very large repulsive forces occur. Choosing a Langevin thermostat 47|, |48!, |49|, the 
system then is equilibrated at T* = 1.16 and is replicated three times in x and y-directions, 
to obtain a system with linear dimensions L^ = Ly = 60. This 9 times larger system then 
is equilibrated again, for a time of 300 MD time units (1.5 x 10^ MD steps), to remove 
the effects due to original periodicity at L^. = Ly = 20. It was carefully tested that for 
the chosen conditions (i.e., for a supercritical solution of relatively short chains) such a 
short re-equilibration time actually was enough. Then the thermostat was switched off, 
and a Galilei transformation of particle velocities was applied to remove the motion of the 
center of mass of our model system. This still rather small system, as described above, was 
only used for testing our simulation and analysis procedures as well as for choosing optimal 
parameters for the pressure-jump simulations. To obtain the initial state of the full system 
at the desired dimensions L^ = Ly = 240, Lz = 12, the system was replicated again 4 
times in x- and y-direction, and the procedure of equilibration and Galilei transformation, 
as described above, was repeated again. The structure factor of the system was carefully 
analyzed to check that any signs of the Bragg peaks (due to the periodic arrangement of 
the replicas) have disappeared. Then the system was equilibrated further for 400 MD time 
units (2 X 10^ MD steps), before the quench was started. Note that at this stage we have 
already a total number of A^ = 589999 particles, namely 294000 chain segments (i.e., 58800 
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chains), 88400 solvent particles and 207599 wall particles. Since each CO2 solvent particle 
contains 3 atoms, and each C16H34 chain contains 50 atoms, the total number of atoms (if 
we had an atomistic model) in our system would be 3205200 (not counting the wall atoms). 

The pressure-jump quenching experiment has the effect that the system after the sudden 
quench can take a larger volume, and since the particle numbers always are fixed, this 
corresponds to a decrease of density. We do not attempt to precisely mimic how the pressure 
jump is carried out in an actual experiment, but we simply rescale the positions of the centers 
of mass of hexadecane and carbon dioxide molecules in three directions such that the final 
dimensions of the simulation box were Lx=Ly=300 and L^ = 15. Of course, one must 
not simply rescale all the coordinates of the effective monomers, since the conformation of 
an individual C16H34 molecule (bond lengths and positions of the monomers along a chain 
relative to each other) should not be rescaled but rather stays the same, just the molecules 
are moved farther apart form each other at lower density. Note that due to Eq. ([3]) this final 
size of the box corresponds to L^ = Ly = 1356 A while L^ = 67.8 A, so the system still is a 
ultrathin nanoscopic film. Wall particles were removed from the system before the rescaling 
of CO2 and C16H34 positions and inserted just after the rescaling procedure, so that the 
arrangement of the wall particles (and the distances between them) stay exactly the same 
as before the quench. In addition, we reduce the energy of the system by rescaling the 
kinetic energy of the particles to ensure that the temperature of the system after the quench 
becomes very similar to the initial temperature of the equilibrated homogeneous system. 
Physically the walls confining a thin fluid film would be massive solid walls of a suitable 
device, of course, and thermostating the walls would have the effect of maintaining constant 
temperature conditions. Our procedure is meant as a short-cut for such a situation. 

The simulation of the system after the quench is performed in the NVE ensemble for 
4000 MD time units corresponding to two million MD time steps. For the first 200 MD time 
units (10^ MD steps), 30 runs were performed in parallel, storing configurations after every 
10 MD time units. For later times, due to the large computational effort for our system, 
five systems were propagated and configurations were analyzed for every 100 MD time units 
only. These simulations were carried out on the multiprocessor system JUMP of the John 
von Neumann Institute for Computing (NIC) at Jiilich, utilizing 16 processors in parallel, 
and the cluster of the SOFTCOMP EU Network of Excellence, utilizing 4 processors in 
parallel. 



III. SIMULATION RESULTS 

A. Transient segregation between solvent and polymer forming a layered state 

Fig. [2] shows typical snapshot pictures that illustrate the time evolution of the phase 
separation process in the thin film. These snapshots show quasi-two-dimensional slices 
parallel to the wall, the left column being about 3 layers (of total 10) away from the wall, 
the right column being close to the center of the film (layer 5). Already these pictures show 
an interesting interplay of phase separation in the directions perpendicular and parallel to 
the confining walls: in the initial stage, t = 10 (Figs. [2^, f), the system is still laterally 
homogeneous, apart from very strongly localized density fluctuations, but there is a strong 
variation of density across the film: most of the effective monomers are concentrated in the 
center of the film (Fig. [2):). This observation is still true at t = 50 (Figs. [2}3,g), but now 
lateral phase separation has clearly started: in the center of the pore (Fig. [2^), the white 
"holes" mean that CO2 bubbles with a few hexadecane molecules (i.e., a dilute solution 
of chains in supercritical CO2) have formed within the concentrated C16H34/CO2 solution, 
while near the walls (Fig. [2)d) we rather have ramified clusters of C16H34 molecules in a 
C02-rich fiuid. At later times, these structures coarsen (t=100, 200) and, at the same time, 
the difference in density between the center of the thin film and the regions near the walls 
diminishes. For t = 1000 (as well as for later times, that are not shown here) the density 
difference has almost vanished (Figs. [2^, j). What is more important, the regions in the 
(x, y)-plane where the CO2-C16H34 interfaces occur, are identical in the left and the right 
snapshot: We can picture the phase separation in the late stages, where the characteristic 
linear dimensions of the growing CO2 bubbles in the concentrated C16H34/CO2 solution in 
X, y- direction are much larger than the film thickness L^, simply as a quasi- two-dimensional 
arrangement of fiat cylinders of height L^, forming bridges between the two walls. While 
some of the C02-droplets at t = 1000 still deviate strongly from a circular cross-section in 
X, y-directions, actually an inspection of snapshots at still later times, such as t = 4000 (not 
shown here) shows that the droplets in fact do develop towards becoming a circular cross 
section, thus, minimizing the interfacial area (and energy). Ultimately (at time t — > cxd in 
a macroscopically large system in x, y-directions) we expect a population of strictly regular 
cylinders of typical radius R{t) connecting the two walls and with R{t) growing to infinity 



as well. Unfortunately, for t = 4000 in our system the number n{t) of cylinders was only 
n{t = 4000) = 8, implying that the data suffer from very strong finite size effects! In 
fact, studies of coarsening in simple diffusive models have suggested that finite size effects 
become important already when the number n{t) of growing domains becomes distinctly 
smaller than n(t) = 20 |57|, and hence our data for t > 1000 clearly suffer from finite size 
effects, despite the rather large linear dimensions and number of particles in our system. 
Therefore, there is no point in carrying out our MD simulations for the system dimensions 
chosen here after longer times. 

Fig. [3] shows now the laterally averaged total density of particles which we define as 
follows 

pUz) = {N^ + Nf)/V, , (7) 

where N^ is the number of solvent (carbon dioxide) molecules in layer z = Zi with Zi 
located in the middle of the interval Az = 0.15; Nf is the number of effective monomers 
of hexadecane in this layer, and Vi = L^LyAz the associated volume of layer i. We have 
tested that the dependence of the profile ptot(^) on the width Az of this volume slice Vi is 
not important. Our choice was taken to ensure that fluctuations in ptot{z) are small enough 
but no significant information on the inhomogeneity of the system in z-direction is lost. One 
can see that near both walls there is always a region (of thickness ~ 0.82cr) essentially free 
of particles, and then the density both in the initial state and in the final state gradually 
increases to an almost constant density in the inner part of the film, ptot{z/Lz ~ 0.5) ~ 0.8 
before the quench and ptot{z/Lz ~ 0.5) ~ 0.4 during the late stages. However, at early 
times after the quench, most of the particles accumulate in the center of the film, so the 
system initially takes a state which shows a phase separation of the liquid-vapor type in 
the z-direction perpendicular to the walls: vapor layers occur close to the walls, and a fluid 
(where CO2 and C16H34 are still almost homogeneously mixed) occurs in the center of the 
film. However, this vertically separated state is unstable against lateral phase separation, 
in the directions parallel to the walls, as we have already seen from the snapshot pictures 
in Fig. [21 Thus, in a particular distance z from the walls the density ptot (-2) approaches its 
final equilibrium in a non-monotonic fashion as a function of time, e.g., for z/Lz = 0.2 the 
density decays fast from a rather large value (ptot ~ 0.6) to a very small value (ptot ~ 0.07) 
immediately after the quench, and only when the lateral phase separation starts the density 
increases again. Note that the two-stage character of the phase separation process, where 

9 



first a stratified structure forms, with low density in the film center, which then laterally 
decomposes, is not a consequence of studying a binary mixture, but rather a consequence 
of purely repulsive wall-particle interactions. In fact, we have checked for pure CO2 that a 
similar behavior occurs. 

It is also interesting to study the density profiles Ps{z) = N^/Vi and Pp{z) = Nf/Vi of 
the solvent particles and the monomers separately, and to define also a profile c{z) of the 
relative concentration of CO2, 

ciz)=ps{z)/[ps{z)+pp{z)] . (8) 

Fig. m displays the latter profile: one can see that despite the fact that we have chosen the 
same repulsive potential between wall atoms and the monomers or solvent particles, respec- 
tively, nevertheless the relative concentration of CO2 near the walls is strongly enhanced, 
both in the initial state and during the late stages of phase separation. Actually, the curves 
for c{z) in the initial state and in the late stages {t > 500) almost fall on top of each other! 
Only in the early times after the quench (10 < t < 100) we do see a much stronger variation 
of c{z): near the walls almost pure CO2 phase is reached. So the phase separation clearly 
proceeds in two stages: induced by the walls, first in the direction perpendicular to the walls 
a layered structure forms, CO2 and C16H34 get almost completely segregated, with the poly- 
mer film in the center and two CO2-C16H34 interfaces near z/Lz ~ 0.25 or 0.75, respectively. 
However, this state costs far too much (interfacial) energy, and is hence unstable towards 
lateral phase separation. Both in the initially homogeneous state and in the final state with 
the "cylindrical" C02-domains across the thin film (Fig. [2]) we have a strong concentration 
enhancement of CO2 near the walls. This enhancement clearly is an entropic effect, from 
the point of view of configurational entropy polymers tend to avoid the regions close to the 
walls. 

Of course, for times t < 100 the system clearly is rather far from equilibrium, and its state 



is changing rather rapidly. Monitoring the ternperature (from the kinetic energy) 47 



48 



and the pressure (using the virial theorem 47|, |48|, |49|]) as a function of time after the quench. 



49| 



a distinct but relatively small increase with time in the region 10 < t < 100 is indeed found 
(Fig. [5]). However, for later times both quantities settle down at constant values, as desired. 
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B. Equal-time structure factors 

In experimental studies of phase separation kinetics, most often the equal-time structure 



factor S(q, t) i s monitored, q being the wavevector of a scattering experiment 2^, 



23 



24 



25 



26 



27 



28l | . For a thin film geometry, q needs to be oriented in the {x, ?/)-plane, of 
course, q = q\\. Also, due to the inhomogeneity of the system in the z-direction (Figs. [3l Hj), 
it is of interest to distinguish in the structure factor from which slice {z) the scattering 
particles contribute to the scattering intensity. Moreover, having two components (which 
here we symbolically denote as A and B, in order to make contact between our notation 



and the relevant literature 



58 



59l |) one must distinguish partial structure factors and those 
which monitor density and concentration fluctuations. Thus, we define the partial structure 
factors, resolved with respect to the 2;-coordinate, as follows 

Sap{q\\,z,t) = -j;^'^'^{ewW- irM{t))]), (9) 

fc=i e=i 

where a,(3 = A or B, fki{t) = fk{t) — r£{t), and N^, Nf^ being the numbers of particles of 

type A or B in the slice at z (i.e., the coordinates Zk{t), zi{t) of the particles must be in the 

range z — Az < Zfc(t), zi{t) < z + Az). While ideally one would like to consider the limit 

/S.Z — > 0, in practice we had to choose a rather larger value of A2; (namely Az = 0.75) in 

order to get enough statistics. 

For fluids Sai3{q\\, z, t) depends only on the magnitude q\\ of q\\ and not its direction. Thus, 

the structure factors monitoring fluctuations of number density (S'nn) and of concentration 

{Sec) are defined as follow, 

Snn{q\\,Z,t) = SAA{q\\,Z,t) +2SAB{q\\,Z,t) + SBBiq\\,Z,t) , (10) 

Scc{q\\,Z,t) = XBSAA{q\\,Z,t) +x\SBB{q\\,Z,t) -2xAXBSAB{q\\,Z,t) (11) 

where xa = Na/{Na + Nb) and xb = Nb/{Na + Nb) are the relative concentrations of 
A(B) particles in the slice centered at z. To simplify the notation in the following, the index 
II from g|| will be omitted. 

Note that due to the motion of particles in the 2;-direction, A^^ and Nb are not conserved 
for a selected layer; in particular, during the early stages of phase separation, xa and xb 
change strongly with time. As an example. Fig. [6] presents S'nn(5', -2, t) as a function of 
wavenumber q for different choices of z and three times: before the quench (a), at t = 40 
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(b) and t = 100 (c) after the quench. The values of z shown in the figure are symmetric 
around the center of the film, which occurs at {Lz/2), and therefore pairs of curves should 
superpose, apart from statistical errors. We see that this symmetry indeed is rather well 
satisfied (e.g., the curves for the layers 5 and 6 are indistinguishable from each other in the 
scale shown in Fig. Et), and hence the statistical errors of our data indeed are rather well 
under control. One can see a peak near q = 2tt which changes relatively little with time: 
this peak and the structure at still larger q reflect the local packing of particles in a dense 
fluid. Apart from the values of z very close to the walls (e.g., for layers 1 and 10, where 
almost no particles occur, as the density profiles p{z,t) in Fig. [3] show), all curves exhibit 
a minimum somewhere in the region 2 < g < 4, while for smaller q, the structure factor 
SnniQ, z,t) increases again. In equilibrium, the maximum of the structure factor occurs for 
g — > (Fig. [6^), as expected, while after the quench for large enough times, Snn{(l,z,t) 
exhibits a well-defined maximum for small q (Fig. ^): this small-angle scattering is the 
"hall mark" of spinodal decomposition. However, close to the walls (i.e., for layers 2, 3, 8 
and 9 centered at 2; = 2.25, 3.75, 11.25 and z = 12.75, respectively, in Fig. [6b) the scattering 
intensity for small q does not seem to decrease again, and so the maximum is much less 
pronounced. Indeed, this range of z clearly exhibits a lack of conservation of the density, 
due to the rapid change of the total density p{z, t) in this regime of times (cf. Fig. [3]). 

The analysis of Scdq, z, t) gives a similar picture, and hence is not shown here. We rather 
try to use both Snn{(l,z,t) and Scc{q,z,t) to extract characteristic lengths R{t) by taking 

suitable ratios of moments ^]. We define Ri{z,t) from Snn{(l,z,t) 

g=gcut q=qcut 

Ri{z,t) = 2nJ2Snn{q,z,t)/Y,(lSnn{q,z,t) , (12) 

q=0 q=0 

and similarly for R2{z, t) from Scc{q, z, t). The resulting data are shown in Figs. [7] and [H] for 
different choices of z and different values for the cutoff gcut- 

Data for 1 < t < 10 were not included, since at such extremely short times after the 
quench both pressure and temperature still are rather strongly time- dependent, the system 
is very far from equilibrium in all respects, and a discussion of the evolution of the system in 



terms of the concepts on coarsening 2^, l23|, |2J, |25|, l26|, |27|, |28j would be rather misleading. 



Also the behavior in the next decade, 10 < t < 100, is difficult to interpret: we see an 
unusually strong dependence of both Ri{z,t) and R2{z,t) on the cutoff gcut, and for some 
values of z there occurs a slight maximum at about 30 < t < 60. This behavior can be 
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attributed to the special interplay between phase separation in the directions parallel and 
perpendicular to the walls (Figs. W^- Since the redistribution of both density (Fig. [3]) and 
relative concentration (Fig. H]) between different z is so pronounced during this range of 
times, the time evolution for a given value of z is similar to the time evolution of a system 
whose order parameter is not conserved. In the thin film as a whole, however, both particle 
numbers are conserved; hence, the density and concentration are conserved variables, when 
we consider the total film. Only for times t = 500 or larger the profiles of density p{z) 
and concentration c{z) are practically independent of time, and then in a particular layer 
(i.e., particular value of z) the order parameters behave as if they were strictly conserved. 
Gratifyingly, in the time region from 200< t < 2000 the data for Ri{z,t) and R2{z,t) show 
indeed a much more standard behavior, being essentially independent from the cutoff gcut, 
and almost independent of z, showing that now a well-defined unique length scale exists in 
the system. Figs. [TJ [8] reveal that in this range of times we almost find straight lines at 
the log- log plots, with a slope slightly below 1/3. For t > 2000 the curves even get slightly 
flatter, so the growth gets slower; we attribute this effect to the onset of finite size effects. 
An important finding of our study, however, is that we do not see any evidence for the 



anomalous law i{t) oc t"^'^ found by Das et al. 29|, ISOj in a symmetric mixture confined in 
thin film geometry. It remains to be understood whether this different coarsening behavior 
is primarily due to the lack of symmetry between phase separating species in our system or 
due to different boundary conditions at the walls. 

C. Equal-time correlation functions in real space 

In the context of simulations, it has some practical advantages to extract characteristic 
lengths from the equal-time correlation functions in real space ^] rather than using the 



structure factors. In fact, also t 
from such a real-space analysis 



le result for a length i{t) growing as i{t) oc t^^^ was extracted 
29l . |30| . Thus, it is of interest to study real-space correlations 
in the present context, too, to check whether the findings of Figs. [7] and [8] are corroborated. 
We define a normalized pair correlation function G{r,z,t) as 

where g{r, z,t) in the equal-time radial distribution function for effective monomers of hex- 
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adecane in a slice z of width 2I\z = 1.5a at the time t. Here, the distance r = \rk{t) — ri{t)\ 
and the coordinates Zk{t), ze{t) of the particles labeled as k and i are restricted to this slice, 
as in Eq. i^. The value of g{r = 0, z, t), which is used to normalize G{r, z, t), we obtain by 
extrapolating g{r, z, t) from the region r > 4cr to r = as described below, thus ignoring the 
local packing effects, which are present in the radial distribution function at short distances. 
Fig. [9] shows such data for the layer 3 {z = 3.75cr) (a) and the layer 5 {z = 6.75(t) (b). 



While in the center of the film 
follow the traditional method 



layer 5) the curves intersect the abscissa, and hence one could 
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30|] to define a characteristic domain linear dimension £{t) 
from the first zero crossing of G{r, z, t), this method clearly does not work in slices close to 
the walls: e.g., for z = 3.75a and time t = 40 (shown in the inset of Fig. ^j^) we rather see 
a continuous decay towards the abscissa with a very fiat minimum at r ^ 18cr instead of a 
clearly identifiable crossing of the abscissa. Thus, we tried heuristically an alternative way to 
extract a length i{t), by fitting a straight line to G{r, z, t) in the regime 0.6 < G{r, z, t) < 1. 
The zero-crossing of these straight lines would allow to identify a length i{t) for all values 
of z. However, this method also is doubtful, particularly for times t < 100, since there the 
curves for G{r, z, t) show strong oscillations for small r. These oscillations are not due to bad 
statistics, but simply reflect the liquid short range order: oscillatory variation of g{r,z,t) 
due to the packing of particles in the nearest neighbor shell, next nearest neighbor shell, 
third nearest neighbor shell, etc., around a particle 58|. This short range order needs to be 
disentangled from the growth of a length scale due to phase separation (Fig. [TOl) . Only when 
g{r) is distinctly nonzero for r > 4cr, the growing length scale can be identified; therefore, 
very short times (such as t = 20) obviously must be discarded. However, in the regime 
r > Aa and t < 100, g{r,z,t) exhibits also clear curvature, and hence any straight line fit 
prone to large systematic errors. Therefore, we choose the ad hoc form 

y{r, z,t) = l + a{z, t) exp[-r/[(2;, t)] + h{z, t) (14) 

to smooth our data for g{r,z,t) using a{z,t), l{z,t), and b{z,t) as fit parameters; b{z,t) 
is negative when a zero crossing (i.e., g{r,z,t) — 1 = at r > Aa) occurs and is positive 
otherwise. This fit-function is used in the range r > 4a and g{r,z,t) > 1.01 but actually 
provides a good representation of the actual data down to and below the first-zero crossing 
for those values of z and t where such a zero crossing occurs. Also, this function is used 
to extrapolated from the region r > 4cr to r = to obtain the normalization constant 
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g{r = 0,z,t) = 1 + a{z,t) + b{z,t) for Eq. (fT3l) . From the first-order Taylor expansion of 
Eq. flT4|) we obtain an approximation to G{r, z, t) 

G(r...t)^ °'--''i\-:f|^:;''-;'''''-'' = i-r/i(..t) (15) 

a{z,t) + o(z, t) 

from which we can also define a characteristic length l{z,t) = l{z,t){a{z,t) + b{z,t))/a{z,t) 
as a value of r at the intersection of a line given by Eq. ( IT5l) with the axis G = 0. 

Fig. [TT] shows the time evolution of the characteristic length i{z, t) extracted in this way 
for three choices of z. While for t < 100, where i{z, t) is only of the order of a few a, indeed 
the data seem to be compatible with a behavior £(z, t) oc t^/^ as observed by Das et al. 29|, 
for the decade 100 <t< 1000 the data seem to be compatible with i{z,t) oc t^/^ showing 
a crossover to i{z,t) oc t^^'^ at later time. For t > 1000, a crossover to a still slower growth 
is evident in our simulations, which is however strongly affected by finite size effects, since 
the number of growing (cylinder-shaped) domains is already rather small. Note that this 
implies that finite size effects already set in for i{z,t) ^ 20a, so despite our large system 
(300cr X 300cr) we cannot follow the kinetics of spinodal decomposition for a large enough 
range of times in order to make significant statements on the asymptotic power law growth! 
Much larger systems need to be simulated for this purpose. 

IV. CONCLUSIONS 

In this work, we have presented computer simulations of spinodal decomposition of a 
coarse-grained model for a compressible binary fluid mixture, which roughly describes hex- 
adecane dissolved in supercritical carbon dioxide. This system has a very asymmetric phase 
diagram in the plane of variables pressure and molar fraction (Fig. [1]), and the pressure- 
jump considered in the present work is strongly off-critical: the number of CO2 molecules 
is 88400 while the number of C16H34 chains is 58800 (leading to 294000 effective segments). 
We find that the phase separation is a two-step process: in the first step, there is a strong 
segregation between solvent and polymer leading to a layered structure, with solvent rich 
layers adjacent to the walls, and a polymer-rich ultrathin polymer film "sandwiched" in be- 
tween. This stratified structure, however, is unstable: the free-standing polymer film in the 
center of the slit pore breaks up, C02-rich bubbles form, and finally a pattern develops with 
cylinder-shaped C02-rich domains, the radius of which grows with a i{t) oc t^'^ law (at the 
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latest stages accessible to our simulation, as long as finite size effects are still negligible). For 
the earlier stages of phase separation, where a strong coupling between the phase separation 
in perpendicular and parallel directions (with respect to the walls) occurs, we conclude that 
a description of the structure in terms of power laws of characteristic linear dimensions is 
somewhat misleading, since characteristic lengths extracted from the structure factor and 
from the pair correlation function are quite incompatible with each other. We suggest that 
there is no simple scaling behavior in this regime. Note that although we have chosen the 
potential between wall particles and solvent particles identical to the potential between wall 
particles and effective monomers, both in the initial and final stages of phase separation 
there is significant enrichment of CO2 near the walls, although our snapshot pictures (re- 
solved as function of z and t in Fig. [2]) indicated that the data still belong to an incomplete 
wetting regime. 

Thus, it would be interesting to have a more detailed theoretical understanding from 
analytical theory for phase separation with two coupled order parameters (density and 
concentration, in our case). Additional simulations would be valuable where wall-particle 
interactions are chosen such that a strictly "neutral wall" situation is achieved, where 
no surface enrichment occurs, and hence the pure confinement effect on phase separation 
(not disturbed by the formation of precursors of wetting layers) could be studied. Also, it 
would be clearly worthwhile to study more systematically how the phase separation kinetics 
depends on slit thickness, composition of the mixture, and quench depth. We did some 
preliminary runs at one different quench depth in which the system volume was increased 
by factor 1.73 instead of 1.95 discussed in our paper and found a rather similar behavior. 
However, significantly different behavior is expected for shallow quenches through the 
critical point, since the correlation length of density and concentration fluctuations can 
exceed the slit thickness, and formation of a stratifled structure is not expected. Finally, 
in order to get rid of finite size effects, simulations with billions of particles on massively 
parallel supercomputers would be required. However, such detailed and computationaly 
extensive studies would be a very challenging task for presently available computer resources 
and itself could be a topic for forecoming papers. We also hope the present work will 
stimulate more work on the kinetics of phase separation in nanoscopic confinement such as 
very thin channels. 
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FIGURE CAPTION 

FIG.l. Isothermal slice through the binary phase diagram of the present model for CO2 
+ C16H34 mixtures at T = 486. 2K (reduced temperature T* = 1.16) in the bulk, using the 
pressure p and the molar fraction x of CO2 as variables. The coexistence curve encloses 
a two-phase coexistence region containing a polymer-rich phase (left) and a supercritical 
CO2 vapor (near x = 1, right). The simulations of quenching experiments discussed in the 



present paper are done for x = 0.6. This phase diagram is taken from the results of Ref. 35|. 



FIG. 2. Snapshot picture showing the structure formation after the quench for 
L^ X Ly = 300a shces of width 1.5cr centered at z = 3.75cr (a-e) and ai z = 6.75a (f-j). 
Snapshots are presented at times t = 10 (a, f), 50 (b, g), 100 (c, h), 200 (d, i) and 1000 
(e, j). The insets in snapshots (c) and (h) illustrate the enlarged regions in the left-bottom 
corner of size 30cr x 30cr (marked by rectangles): the gray spheres correspond to the 
supercritical solvent molecules, and the black ones represent the chain molecules. 

FIG. 3. Total density profile ptot(^) (laterally averaged in the thin film) plotted vs. 
position across the slit pore z/Lz for the initial state before the quench (L^ = 12cr, 
Lx = Ly = 240cr), and for different stages of phase separation (L^ = 15cr, L^ = Ly = 300cr). 
The curves are shown at different times after the quench, as indicated. 

FIG. 4. Relative concentration profile c{z) of the solvent laterally averaged in the thin 
film plotted vs. position across the slit pore z/L^ for the same times and the system 
dimensions as shown in Fig. [3l 

FIG. 5. Time evolution of the reduced temperature T* (top) and the reduced pressure p* 
(bottom) before and after the quench. The time of the quench is at t = 0. 

FIG. 6. Density-density structure factor Sji^{q,z,t) in the initial state after equilibration 
(just before the quench) for the system of the linear dimensions L^ = Ly = 240cr, Lz = 12a 
(a), and after the quench at times t = 40 (b) and t = 100 (c) for the system dimensions 
Lx = Ly = 300(7, Lz = 15cr. Various choices of z are included, as indicated in the figure. In 
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case (a), an average over 130 configurations was performed, wliile in cases (b) and (c), 30 
configurations were averaged over. Note: tlie two curves for layers 5 and 6 (in tlie middle 
of tlie slit pore) cannot be distinguished on the scale shown in figure (c). In all cases, the 
film was divided into 10 slices of width 2Az = L^/IO. 

FIG. 7. Log-log plot of the characteristic domain size Ri{z,t) (calculated from the 
density-density structure factor ^nn shown in Fig. [6]) vs. time for four choices of z, as 
indicated in the figure, and three plausible choices of the cutoff for the wave vector: 
gcut = 1,2, or 3, respectively. 

FIG. 8. Log-log plot of the characteristic domain size R2{z,t) (calculated from the 
concentration-concentration structure factor Sec) vs. time for four choices of z, as indicated 
in the figure, and three plausible choices of the cutoff for the wave vector: g^ut = 1,2, or 3, 
respectively. 

FIG. 9. The real-space normalized pair correlation function G{r,z,t) [Eq. ( TT3l) ] vs. 
distance r for z = 3.75a (a) and z = 6.75a (b) at various times, as indicated in the figure. 
The inset in part (a) shows a magnification of the region near the origin for early stage of 
spinodal decomposition. Note, that the curve for t = 40 exhibits a very fiat local minimum 
near r ^ 18cr, which prohibits using the first-zero crossing of G{r,z,t) to measure the 
characteristic domain size in this range of the time. Straight lines illustrate fits to estimate 
the characteristic domain length scale l{z,t) from an effective initial slope of these curves 
at r = using Eq. (fT5|) . 



FIG. 10. Plot of the radial distribution function g{r, z, t) for the middle layer {z = 6.75a) 



vs. distance r for various times, as indicated. 



FIG.ll. Log-log plot of characteristic domain linear dimension vs. time calculated in 
three layers of thickness Az = 1.5a parallel to the walls, centered at 2; = 3.75cr, 5.25a and 
z = 6.75a, respectively. The domain size shown here is defined in the text by Eq. (fT5l) . 
Straight lines are the guide for eyes illustrating the power laws £{t) oc t^^'^ (dashed line), 
i{t) oc t^/^ (dashed-dotted line) and i{t) oc t^^^ (dashed-double dotted line), respectively. 
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